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Section  1 


The  Method  of  Frieman  and  Kroll 


The  derivation  of  the  method  developed  by  Frieman  and  Kroll  (1973)  for  the  calcula¬ 
tion  of  electromagnetic  fields  due  to  a  transmitting  antenna  within  a  one-dimensional  profile  of 
electrical  conductivity  starts  with  Maxwell’s  equations  in  the  charge-rationalized  MKS  system  of 

units; 

VxW- 7+ 

ot 

(1.1) 

Vx£--M 

dt 

(1.2) 

VD  —  p 

(1.3) 

v-5-o 

(1.4) 

Assuming  a  time  dependence  of  eml  for  the  fields  and  currents,  the  first  two  equations  become 

Vx£  —  —iuifjiH 

(1.5) 

Vx#-  J  +  iuieE  , 

(1.6) 

where  the  linearity  conditions 

D  =  *E 

(1.7) 

B-nH 

(1.8) 

are  assumed,  with  the  electrical  permittivity  «  and  the  magnetic  permeability  fi  taken  to  be 
independent  of  the  field  variables  and  not  explicitly  dependent  on  time.  The  choice  of  e"“'  over 
e~'"“  as  the  time  dependent  factor  is  motivated  by  the  fact  that  this  choice  makes  it  possible  to 
directly  compare  the  field  values  calculated  by  this  method  with  the  Fourier  coefficients  of  field 
value  measurement  real  time  series  that  have  been  Fourier  analyzed  according  to  the  standard 
convention.  If  one  then  expresses  the  current  density  distribution  7  in  the  form 


J  ”  J4  +  rrE 


(1.9) 
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where  JA  is  the  current  uensity  distribution  in  the  transmitting  antenna  and  the  other  term  on 
the  right  is  the  current  density  distribution  induced  by  ohmic  processes  in  material  with  an 
electrical  conductivity  distribution  it,  then  (1.6)  becomes 


Vx//  =  JA  +  E 

ItUfJL 

y 2  “  iiundiot+a) 

(1.10) 

(l.ii) 

where  m.  «.  and  or  are  functions  only  of  position  in  space. 

A  vector  potential  A  is  now  defined  such  that 

//  =  Vxl 

(1.12) 

E  =  i<ufi  —  Al 

y  \ 

(1.13) 

If  n  is  taken  to  be  a  constant,  then  equation  (1.5)  is  satisfied  trivially;  accordingly,  henceforth 

fi  is  taken  to  be  uniformly  equal  to  its  free  space  value  /xo-  Working  with  equation  (1.10)  gives 

Vx//=  Vx(VxA) 

(1.14) 

-  V(V-J)  -  V2?  *  J4  +  y2|v-yVJ  _ 

V2I  -  y'A  -  V(V  A)  -  y2V^V-J  =  -1A 

(1.15) 

y2v\v-A  -  V(V-I)  -  (V  ^)-Vv(r) 

y  r 

(1.16) 

V}A-y2A-  (V-.4)~V(y2)  = -J4  . 

y 

(1.17) 

Equations  (1.12)  and  (1.13)  admit  a  restricted  gauge  transformation  for  A,  such  that 
A'  —  A+V>li  and  A  generate  the  same  fields  if  V<b  satisfies  equation  (1.17)  with  JA  set  uni¬ 
formly  equal  to  zero. 

If  «  and  <r  are  taken  to  be  functions  only  of  the  altitude  coordinate  z.  then  equation 
(1.17)  becomes 

V2^  -  y2A  -  :(V'I)i^  —7,  .  (1.18) 

y*  CiZ 

whose  homogeneous  solutions  can  be  divided  into  solutions  with  A  in  the  r  direction  and 
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solutions  with  A  perpendicular  to  the  z  direction.  In  the  first  case  we  assume  the  solutions  to 
be  separated  in  the  form 

A  -  f(z)F(x,y)z  ,  U.19) 

which  gives  from  (1.18)  with  set  equal  to  zero 

-  -V  &  4£(?--  -  (y2+k2)f(z)  °  0  (1.20) 

dz2  y2  dz  dz 

V2HF(x,y)  +  k2F(x.y)  =  0  ,  <1.211 


where  the  subscript  H  indicates  that  derivatives  occur  only  in  the  horizontal  plane.  From  sym¬ 
metry  considerations  one  has  that  the  solutions  arising  from  the  separation  (1.19)  can  give  only 
transverse  magnetic  fields.  Boundary  conditions  on  f (z)  may  be  deduced  from  the  facts  that  the 
transverse  field  components  Er  and  Hr  must  be  continuous  across  discontinuities  in  o\  and 
£r— 0  at  a  discontinuity  on  one  side  of  which  or=*=°o.  From  (1.13)  we  have  for  £r 

Et  -  icjn\^^-VHF(x,y)  .  (1.22) 

UZ 


which  gives  that 

1  dfx(z) 1  df2(z) 

y 2  dZ  yl  dz 


(1.23) 


across  a  discontinuity  in  o-(z)  at  z— z0,  and 


df(z) 
dz  :° 


0 


(1.24) 


at  a  discontinuity  at  z-z0  in  cr(r)  on  one  side  of  which  o-=*°o,  and  from  (1.12)  we  have 


Ht  —  /(z){Vx£(x.y)z) 


(1.25) 


from  which  we  have 


f\(z0)  -/2(z0)  (1.26) 

across  a  discontinuity  at  z—zo  in  o-(z).  In  the  case  of  horizontal  homogeneous  solutions  for  A 
we  assume  separation  in  the  form 

A  -  <j(z)Q(x,y)  .  (1.27) 

where  (5  has  no  z  component  and  we  exploit  the  previously  mentioned  gauge  freedom  to  make 
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the  constraint 

V(5  =  0  (1.28) 

This  separation  gives 

-  (y2+k2)q(z)  -  0  (1.29) 

az 

VhQ{x,y)  +  mx,y)  =  0  (1.30) 

Note  that  by  (1.27)  and  (1.28)  we  have 

VA  =  q(z)\7Q  +  QVq(z)  =  0  ,  (1.31) 


so  according  to  (1.13)  this  class  of  solutions  to  (1.18)  can  give  only  transverse  electric  fields. 
For  Ej  we  have  from  (1.13)  and  (1.27) 

Et  —  -i(t>nq(z)Q(x,y)  ,  (1.32) 

which  implies  that 

q\(za)  -  qiiz0)  (1.33) 

across  a  discontinuity  at  z=zo  in  cr(z),  and 


q  (r0)  =  0 

at  a  discontinuity  at  z«=zo  in  <r(z)  on  one  side  of  which  cr=~°°.  For  Hr  we  have 
HT  =  ^p-{zxQ{x,y))  , 

which  gives  that  across  a  discontinuity  in  <r(z)  at  z=zo 

dq\(z)  dq2(z) 
dz  :°  dz  :° 


(1.34) 


(1.35) 


(1.36) 


Mutual  orthogonality  of  different  solutions  to  (1.20)  under  identical  boundary  condi¬ 
tions  is  demonstrated  as  follows; 


d_ 

dz 


fm\fn 

y 


y  y 


LhLf  +  f  -Lr  > 

4  fa  J ' 1  '  n 


(1.37) 


± 

dz 


y  y 


(1.38) 
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1  dy' 
T  dz 


=:~(yl+k~)f„ - ^(y-+k 

r  r 

J mJ n  t i  ,  i  v 
"S  '  Aw  km  I 

y‘ 

Ui  ,  }\  C  f nxf n  .  {J mf n  ./ nJ n i  )  I.. 

„‘-A4>  J  — 2  *  - i - /  =  o 

I  y  r  1 


(1.39) 


indicating  that  if  the  eigenvalues  of  the  two  solutions  are  different  then  the  integral  must  be 
equal  to  zero.  The  corresponding  demonstration  for  solutions  to  (1.29)  is 


d_ 

dz 


(dm  Qn  Qn  dm  ^ 


Qm  Qn  Qn  Qm 


(1.40) 


“  dmQnikj-k,;) 


II 

(k2-k2)  J  qmq„  dz  «  (qmQ„'-q„Qm')  /'  =  0 


(1.41) 


If  m— n,  it  is  clear  that  the  integrals  will  always  be  positive  provided  that  the  functions  within 
the  integrals  are  not  uniformly  zero. 

In  anticipation  of  a  future  need.  Green's  functions  solving  the  equations 


« -  y(!),  dz  8; 


-  {y(z)2+A-2},?t  (z,z'.k)  =*  8(r-c') 


(1.42) 


and 


—j gH(z,z\k)  -  {y(z)2+k2)gH(z.z\k)  =  S(z-c’) 


will  now  be  constructed.  A  valid  solution  to  the  first  of  these  equations  is 
,  ,,,  7tO:<,A)/t,(r>,A) 

Wv{k)  -  Ji  (r,A)/;/(r,A)  -  fL'(z,k)fL  (z.k)  , 


(1.43) 


(1.44) 

(1.45) 


and  a  valid  solution  to  the  second  is 
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8h(: 


A)  = 


qL  (z<,k)qL  (tx.A ) 
W„(k) 


(1.46) 


WH(k)  =*  qL(z,k)qL'(z,k)  -  qL'(z.k)qL  (z.k)  ,  (1.47) 

where  /t  and  /c  are  solutions  to  (1.20)  with  A2  specified  and  their  first  derivatives  constrained 
to  vanish  at  the  lower  and  upper  ends  respectively  of  the  range  of  z.  qL  and  </<  are  solutions  to 
(1.29)  with  A2  specified  and  their  values  at  the  lower  and  upper  ends  respectively  of  the  range 
of  z  constrained  to  vanish,  r <  is  whichever  of  z  and  is  the  smaller,  and  is  whichever  of  z 
and  is  the  larger.  Note  that,  if  k?„  is  one  of  the  eigenvalues  of  (1.20)  with  all  constraints 
enforced  and  and  conditions  imposed  to  require  discreteness  of  the  eigenvalues  and  A/J„  is  one 
of  the  eigenvalues  of  (1.29)  under  the  same  circumstances,  g,  is  singular  wherever  k  is  equal 
to  any  of  the  kin  and  gn  is  singular  wherever  k  is  equal  to  any  of  the  kH„.  These  solutions  to 
(1.42)  and  (1.43)  are  useful  for  practical  calculation  of  and  glt  respectively,  but  inappropri¬ 
ate  for  the  purposes  of  this  derivation. 

The  derivation  of  the  alternate  form  of  gt  starts  with  the  identity 


d_ 

dz 


fn\sr  ~ 

r  y* 


^-(A-’-A-r,,) 

r 


(1.48) 


derived  in  a  manner  similar  to  that  which  produced  (1.38),  where  /„  is  /„(;),  g i  is  (c.r'.A), 
y2  is  y2(c),  and  primes  denote  differentiation  with  respect  to  z.  Integration  of  (1.48)  with 
respect  to  z  gives 


0  =  (A2— Af„)  / 


/n(z)gt  (z.z'.k) 
yHz) 


dz  + 


y  :(=’> 


(1.49) 


from  which  we  have 


gy(:.:\k) 

U 

A,  “  f 


1 


( A2„— A2) 
r  fn2(z) 


y2(z ) 


dz 


fn(z')  fn(z) 
y2(z')  /,„ 


(1.50) 

(1.51) 


where  C  is  anything  that  is  orthogonal  to  fn(z)  according  to  the  definition  given  by  ( 1.39);  a 
logical  guess  as  to  the  nature  of  C  gives 


giiz.z'.k) 


1 


(kL~k2) 


y2iz')  />,„ 


(1.52) 
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since  from  (1.39)  we  have 


J,  y:(r) 


d-  S  m„  J  j  ft 


Plugging  (1.52)  back  into  (1.42)  gives 


h  rV)/,„ 


') 


(1.53) 


(1.54) 


exactly  the  expression  produced  if  one  does  a  series  expansion  of  8(r-r’)  in  the  functions 
/„(-).  The  derivation  of  the  alternate  form  of  gH  follows  the  same  pattern,  with  the  result 


=»  1 

k)  =  T  - I - 

(kf,n-k2) 


q,„(:')q,„(z) 


I 


Hm 


(1.55) 


U 

lHn,  =  /  dz 

I 


(1.56) 


Note  that  the  summation  in  (1.52)  starts  with  m=0.  whereas  the  summation  in  (1.55)  starts 
with  m  =  l;  due  to  the  different  boundary  conditions  on  the  solutions  of  (1.20)  and  (1.29)  there 
is  a  0th  order  solution  to  (1.20)  whereas  the  lowest  order  solution  to  (1.29)  is  lsr  order.  Pro¬ 
vided  that  the  /„(:)  and  the  ?n(r)  each  form  complete  sets,  which  they  should  under  the 
proper  boundary  conditions,  expressions  (1.52)  and  (1.55)  should  be  valid  solutions  to  (1.42) 
and  (1.43)  respectively;  one  condition  that  should  give  complete  solution  sets  is  that  <r(r)  goes 
to  °°  at  the  upper  and  lower  ends  of  the  range  of  z. 

In  considering  the  electromagnetic  fields  resulting  from  mode  excitation  by  a  com¬ 
pletely  vertical  current  distribution,  we  define 

J~  =  zJ:  .  (1.57) 


and  from  symmetry  considerations  we  can  write 
.4  =  :A.  . 


(1.58) 


We  define 


•ao 

A.(p.9,:)  -  k  dk  h„,  (z.k)e"""jm(kp) 

m  0 


(1.59) 


its  inverse  operation 


h,„  (z.k)  = 


(1.60) 


8 


x  A_ApJ.z)e~'mHJ,n(kp) 


the  expression 


30 

JAp.O.z)  =  £  f  A  dk  Jm(z,k)em*JJkp\ 

m  *1) 


(1.61) 


and  its  inverse 

,  2n  » 

>-<=.*>  -  4r  /  J 


0  0 


p  dp  d0 


x  JAp,9.;)e~'mHJm(kp) 


(1.62) 


where  Jm(kp )  is  a  Bessel  function  of  the  first  kind  of  order  m.  From  (1.18)  we  have 


V2/!.  —  y2A: 


1  dy~  dA: 
y 2  dz  dz 


(1.63) 


given  the  identity 


V2e‘mHJm(kp )  -  -k2e,mliJJkp) 


(1.64) 


we  substitute  (1.59)  and  (1.61)  into  (1.63)  to  get 

£  J  I  A-  dk  Jm(kp) 


m  0 


&  hm  /  2  ,  i  i  1  c/*y" 

1?-  •  (y  +t')h-  -  7 

Integration  with  over  all  0  gives 


+  Jm 


=  0 


OO 

/ 


A  dk  Jn(kp) 


ft  hn  ,  2 ,  i  f  1  dy‘  ,  ", 

!?•  -  -  7T17  + 


(1.65) 


(1.66) 


for  every  integer  n;  that  (1.66)must  hold  for  all  positive  values  of  p  suggests  that  the  expres¬ 
sion  in  the  brackets  may  be  uniformly  equal  to  zero,  and  indicates  that  even  if  it  need  not  be 
uniformly  equal  to  zero  it  can  safely  be  constrained  to  be  that  way.  as  any  solution  for  h„  that 
satisfies  (1.66)  is  acceptable.  Therefore,  we  can  write 
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h„" - _  h 2+k2)hm  =  -J„  ,  (1.67) 

y*-  az 

where  a  prime  denotes  partial  differentiation  with  respect  to  z,  and  use  of  the  Green's  function 
Si  defined  by  equation  (1.42)  gives  from  this 

14 

tim(z,k)  =■  -  f  gi  (z,:\k)Jm  (z\k)  dz'  (1.68) 

i 

An  alternate  expression  for  A:  is 


A:(p,9,z)  «  X 


30 

Am  (p,z)  -  f  kdk  hm  ( z,k)Jm  ( kp )  ; 


this  last  expression  way  be  rewritten  as 


oo 

A.„ (p,z)  -  [h  J  k  dk  hm(z,k)Ht2i  m(kp) 


(1.69) 

(1.70) 


(1.71) 


by  way  of  the  identity 


J  k  dk  hm(z,k)H^2)  m(kp) 


(1.72) 


—  2  f  k  dk  km(z,k)Jm(kp) 

Jo 

where  Hi2\ (x)  is  a  Hankel  function,  composed  of  Bessel  functions  of  the  first  and  second 
kind  according  to  the  definitions 

H{Um(x)  ~Jm(x)  +  iN„(x)  (1.73) 

ff2)m(x)  «/„(*)  -  iNm(x)  ;  (1.74) 

the  functions  //*2)mUcp)  were  chosen  over  the  functions  J„{kp)  to  describe  the  behavior  of 
the  A:m  with  respect  to  p  because  of  properties  of  their  asymptotic  behavior  whose  utility  will 
become  apparent  below. 

The  identity  (1.72)  derives  from  the  identities 


Jm(-x)  -  (-l)m  Jm(x) 


(1.75) 


HiVm(-x)  -  -(-I)"1  ff2>m(x ) 


(1.76) 
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expression  (1.71)  becomes 


At  \  1  T-  fn(~^  .  .  , 

A;,n  (P'- )  *  j  2^  .  *~nm\fj) 

4  77  „-0  'In 


-  j_ * *  — ,-n: py —  • 

The  Hankel  function  used  in  the  integrand  of  this  integral  has  the  assymptotic  form 

H'!,~{x) -f! 


where  x»  1  and  x»m  , 


(1.84) 

(1.85) 


(1.86) 


which  suggests  that  (1.85)  may  be  solved  by  contour  integration,  with  the  contour  being  drawn 
through  the  lower  half  plane  and  passing  below  the  branch  point  in  Hi2)m(kp)  at  k^O;  the 
solution  is 

L„m(p)- inKn„(kyn)Hi2)m(kvnp)  ,  (1.87) 

where  kv„  is  the  root  of  ky„  that  has  a  negative  imaginary  part.  Hence,  we  have  for  A:m 


A:m(p,z)  =  £  Cnm  fn(z)H{1\(kv„p) 

n—  0 


/  Knm  (kyn) 
4  hn 


(1.88) 

(1.89) 


Equation  (1.86)  indicates  that  the  modes  which  attenuate  most  slowly  with  increasing  distance 
from  the  source  are  those  that  have  the  smallest  values  of  |//m(Av*)I:  accordingly,  the  closer  to 
the  transmitting  antenna  the  point  at  which  one  wishes  to  calculate  the  field  values  is,  the  more 
modes  are  required  to  give  answers  of  acceptable  accuracy.  Note  that  (1.88)  is  a  proper  solution 
to  the  problem  only  if  the  eigenvalues  are  discrete. 

As  an  example,  suppose  that  the  transmitter  is  a  vertical  point  current  dipole  of  dipole 
moment  D,  a  dipole  whose  length  dl  is  made  arbitrarily  small  and  the  current  flow  I  through 
which  is  made  correspondingly  large  in  such  a  way  that  the  dipole  moment  D  of  the  dipole, 
defined  by 


D  -  I  dl 


(1.90) 


has  a  finite  value  in  the  limit  of  infinitely  small  length  and  infinitely  large  current:  then  we  have 
for  J- 
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j.-d  , 

2  irp 

(1.91) 

where  c'  is  the  antenna's  vertical  position  and  its  lateral  coordinate  is  p=0. 

Since  J.  has  no 

9— dependence,  only  terms  of  order  m=0  contribute.  This  being  the  case,  we  i 

from  (1.83)  and  (1.89) 

mmediately  have 

r  iD  Jn(z') 

"°  =  4  y2(z')lyn 

(1.92) 

Cnm  =0  for  m^O 

(1.93) 

In  considering  the  case  of  a  horizontal  current  distribution,  we  take 

J:  —  o  ,  J4  =  JH 

(1.94) 

and  assume  infinitely  conducting  plates  at  the  top  and  bottom  of  the  range  of  the  altitude  coor¬ 
dinate;  in  order  to  maintain  consistency  with  the  behavior  of  the  homogeneous  solutions  at 

such  boundaries,  we  require  that  -4W=0  and  A.'= 0  at  these  upper  and  lower  boundaries,  con¬ 
straints  that  serve  to  fix  the  gauge  and  thereby  permit  a  unique  solution  for  A  everywhere 

between  the  plates.  We  define 

A  =*  Ah  +  zA:  ,  zxAh^O  , 

(1.95) 

which  gives  from  (1.18) 

V2/4 h  -  y2AH  -  -Jh 

(1.96) 

t-t2  ..  1  dy1  dA:  2  1  dy2  „  t 

V/l--  SZ  • 

(1.97) 

We  define  the  expressions  Tm(z,k)  and  Jm(z,k )  such  that 

oo 

IH(p,0,z)  -Zf  k  dk  -sm{z,k)em"J„{kp) 

m  0 

(1.98) 

.  2jt  oo 

f  f  p  dp  dQ  AH(py9,z)e~'m*Jm{kp) 

Jo  Jo 

(1.99) 

oo 

7„ip.0,z)  -  Z  f  k  dk  J„(z,k)e"”HJm(kp) 

m  0 

(1.100) 

J,„(z,k ) 


2^  /  /  PdpdO 

7)  Jo 


(1.101) 
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x  JH(p,0,z)e'm,,Jm(kp)  j  ; 

substitution  of  (1.98)  and  (1.100)  into  (1.96)  with  the  use  of  the  identity  (1.64)  gives 

V'L-./c)  -  (y2+k2)Tm(z,k)  =  -Jm(z,k)  (1.102) 

by  the  same  argument  that  gave  (1.67).  Calculation  of  V-Atl  from  (1.98)  gives 

QO 

VAH  -  £  /  k  dk  ~sm(z,k)-V e'm*Jm (kp)  ,  (1.103) 

m  0 

since  we  have 


VT,(;,()=0  , 

as  Tm  has  no  z  component.  This  motivates  the  expression  of  A:  in  the  form 
AAp,0,z)  =  Z  J  k  dk  pm(z,k)-Vem*  JJkp)  , 


m  0 


2tt  oo 


pm(z.k)  -  — 


2nk2  i, 


p  dp  dO 


(1.104) 


(1.105) 


(1.106) 


x  /4,(p,fl,c)Ve-""HyOT(Arp) 

Substitution  of  (1.103)  and  (1.105)  into  (1.97)  with  use  of  (1.64)  then  gives 

pm"iz,k)  -  \^f-pm'(z,k)  -  ( y2-k2)pjz,k )  (1.107) 

y  uZ 

-  -7 ~t-Jm{z,k) 

y  QZ 

The  requirements  that  -4p“0  and  A:'=*0  at  the  upper  and  lower  boundaries  places  the  same 
conditions  on  7m  and  ~pm  respectively,  giving  7m  the  same  boundary  conditions  with  respect  to  z 
as  qm  and  pm  as  fm. 

Use  of  the  Green's  functions  gH  and  g\  give  for  7m  and  p,„ 


K 

sm(z,k)  -  -  J  gH(z,z\k)  JJz'.k)  dz' 

I 

pm(z,k)  -  Jf  |  gi  (z,z',k) 


(1.108) 


(1.109) 
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Work  with  (1.109)  presents  a  complication,  in  that  both  gt  and  Jm  have  numerous 
singularities  in  the  complex  k  plane;  in  order  for  contour  integration  to  be  practical  these  func¬ 
tions  must  be  expressed  as  series  expansions  about  their  singular  points,  but  if  both  functions 
are  so  expressed  the  integration  process  becomes  impractically  cumbersome.  The  problem  is 
solved  by  defining  two  versions  of  pm ,  in  each  of  which  one  of  the  two  functions  is  expressed 
in  expanded  form  and  the  other  is  left  in  functional  notation; 


PHm^,k)  -  - 


2tt 


Knm(k) 

iHn 


1 

(kf]n-k2) 


un(z,k) 


(1.120) 


u„(z,k)  =■*  J*  gy(z,z',k) 


1 

y2(zj 


dy2(z') 

dz 


(1.121) 


x  q„(z')  dz' 


where  (1.110),  the  series  expansion  of  sm ,  is  used,  and 


Pvm(z,k)  -  - 


Mnm(k) 


lit 


/, 


Vn 


( kl 


Vn' 


- k 2) 


f„(z) 


dy2(z) 
yHz)  dz 


(1.122) 

(1.123) 


x  gH(z,z\k)e  imllJ„(kp)JH(p,9,z) 


where  (1.52),  the  series  expansion  of  is  used  with  (1.108)  and  (1.101).  Although  in  princi¬ 
ple  Pvm  and  pHm  are  each  equivalent  to  the  standard  contour  integration  algorithm  will  not 
see  all  of  the  singularities  of  either,  and  for  purposes  of  contour  integ  ation  they  neatly  divide 
the  singularities  of  pm  between  them;  the  function  pTm, 

PTm  *  Pvm  +  PHm  •  (1.124) 

is  not  equal  to  p„(z,k),  but  within  a  contour  integral  it  will  behave  as  if  it  were.  The  expres¬ 
sions  (1.121)  and  (1.123)  that  contain  the  hidden  singularities  are  simple  enough  in  form  so 
that  contour  integration  will  usually  not  be  required  to  evaluate  them;  expressions  (1.44)  and 
(1.46)  may  be  used  for  the  Green’s  functions  in  the  integrands  of  these  expressions. 

We  can  now  represent  (1.105)  in  the  form 

A.(p,9,z)  -  £  \A:ym(p,9,z)+A.H„(p,9,z)) 


(1.125) 
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/Wm(p,fl,2)  -  f  k  dk  pVm(z,k)Ve‘mHJn,(kp)  (1.126) 

Jo 

oo 

A Mm(p,0,z)  =  f  *  dk  pHm(z.k)Ve"”HJm(kp)  ,  (1.127) 

i) 

with  the  understanding  that  the  integrals  will  be  evaluated  with  contour  integration  .  Inspection 
of  (1.109),  (1.99),  and  (1.52)  shows  that  pm(z,k),  and  hence  pim(z*k)  and  pHm(z,k),  have  the 
same  inversion  properties  with  respect  to  k  as  Jm(kp ),  so  the  same  line  of  argument  that  gave 
(1.114)  and  (1.71)  permits  us  to  write 


x  p„m(z,*)-Vc™*//<2>m(A:p)  j  . 

Writing  these  cut,  we  have 

AzHm(PS,z)  "  -T-T  £  T"  J  kdk  (1.130) 

W  n-l  ‘Hn 

x  -£T(?2 \»n(z,k)-Ve"”»H'»m(kp) 

\kHn~k  ) 

A.ym(p,9,z)  -  ~Y~  £  -J—  J  k  dk  (1.131) 

„_o  hn  -OO 

*  )  • 

with  only  the  poles  explicitly  shown  in  the  denominators  being  evaluated.  The  solutions  are 

4.Wm(p,0,z)  -  £  un(z,kHn)  (1.132) 

rt-1 

X  DnmVe'”",Hi2)JkHnP)  I  , 


where  D„m  is  as  given  in  (1.1 19),  and 


n  -0 


X  Bnm  •  V e‘mUH,2'm  (ki„p) 


(1.133) 


where 


4  A* 


(1.134) 


A  useful  practical  example  for  this  work  is  that  of  a  horizontal  point  current  dipole  of  dipole 
moment  M;  this  is  represented  by  the  expression 

7„(p,0,z)  -  M  8(P)8(z- ;,)  ,  (1.135) 

2np 

where  h  is  a  unit  horizontal  vector  and  z  is  the  altitude  coordinate  of  the  dipole,  which  is  taken 
for  convenience  to  have  the  lateral  coordinate  p-0.  As  with  the  vertical  point  dipole  case,  only 
the  m— 0  order  contributes,  and  the  appropriate  coefficients  are 

D„m  =  --j  Q„(z')h  (1.136) 

4  ‘Hn 


Bn, 

Dm 


B„ 


0  for  m 5*^0 


(1.137) 

(1.138) 


If  one’s  conductivity  profile  cr(z)  is  taken  to  be  composed  of  homogeneous  slabs,  then 
equations  (1.20)  and  (1.29)  have  simple  analytical  solutions  within  each  of  the  slabs,  and  these 
solutions  are  easily  patched  together  at  the  interfaces  by  use  of  the  appropriate  boundary  condi¬ 
tions.  However,  the  solutions  contain  exponential  factors,  and  as  a  practical  consequence  of  this 
the  solution  process  is  numerically  unstable.  Nevertheless,  excellent  results  were  obtained  by 
working  a  Ricatti  transformation  on  the  equations,  and  then  deriving  solutions  to  the 
untransformed  equations  from  the  Ricatti  equation  solutions.  Both  (1.20)  and  (1.29)  can  be 
represented  in  the  form 

^"(x)  +  a|(x)^’(x)  +  floOrb'Or)  "  0  ;  (1.139) 


the  Ricatti  transform  is 
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v(x)  =  I~l  ,  (1.140) 

,v(x) 

and  equation  (1.139),  transformed  by  this  substitution,  becomes  the  Ricatti  equation 

v'(x)  +  v’(.v)  +  a\{x)  v(x)  4-  0O(x)  =  0  •  (1.141) 

From  (1.140)  we  have  that  the  inverse  transform  is 

X 

y(x)  —  C  exp{J*  v(i)dt\  ,  (1.142) 

where  C  is  some  arbitrary  constant;  once  a  solution  to  equation  (1.141)  is  obtained,  the  relation 
(1.142)  can  be  used  to  give  from  it  a  solution  to  the  equation  (1.139).  The  most  general  form 
of  the  Ricatti  equation  is 

v'(x)  +  ai(x)v1(x)  +  0|(x)v(x)  +  a0(x)  =  0  ,  (1.143) 

and  it  is  a  property  of  this  equation  that  if  two  particular  solutions  to  it,  v’i(x)  and  w2(.v) .  are 
known,  then  a  general  solution  v(x)  is  given  by 

«  K  exp{J*a2(r)[v2(/)-v’i(/)]<ff}  ,  (1.144) 

where  K  is  some  arbitrary  constant.  In  the  event  that  the  coefficients  of  (1.143)  are  all  con¬ 
stants,  such  that 

v'(x)  +  av2(x)  +  bv(x)  +  c  —  0  ,  (1.145) 


then  there  are  two  constant  solutions  that  can  be  determined  by  the  quadratic  formula. 


Within  a  homogeneous  slab  both  equations  (1.20)  and  (1.29)  give  the  constant  coefficients 

0—1  ,  b-  0  ,  c — (y2+fc2)  ,  (1.147) 

for  which  the  constant  solutions  are 

vt  —  — /3  ,  v2  -  13  ,  0  —  y/y2+k2  (1.148) 

This  gives  from  (1.144) 


v(x)+0 

v(x)-0 


K  e ^ 


(1.149) 
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Some  algebraic  manipulation  of  this,  with  appropriate  choice  of  K.  gives 


v(x) 


€  =  '/’  lOg 


tanh(j3(x-.v0)+«l 

vo+/3 


Vo— /3 


where 


(1.150) 

(1.151) 


v0=*  v(x0)  ■  (1.152) 

The  corresponding  solution  for  y(x),  developed  from  (1.150)  by  way  of  (1.142),  is 

y(x)  -  y0fcosh{/3(.v-xo))+(v(y'/3)sinh(/3(.v-.vo))l  ,  (1.153) 

where 

v0  -  y  (x0)  .  (1.154) 

The  matching  conditions  at  the  interfaces  between  slabs,  derived  from  (1.23),  (1.26), 
(1.33),  and  (1.36)  by  way  of  (1.140),  are 

(1.155) 

vt  yi 

for  the  Ricatti  solution  to  (1.20),  and 


vi  —  v2  (1.156) 

for  the  Ricatti  solution  to  (1.29).  The  boundary  condition  on  v  at  an  interface  with  an  infinitely 
conducting  plate  on  one  side  is  v— 0  for  the  first  case,  but  for  the  second  case  v  goes  to  infinity 
at  the  boundary,  the  approach  that  was  used  to  deal  with  the  infinitely  conducting  end  plates 
was  to  use  (1.153)  in  combination  with  either  condition  (1.24)  or  condition  (1.34)  to  derive  a 
value  of  v  at  the  side  of  the  end  slab  opposite  to  the  side  interfacing  with  the  infinitely  conduct¬ 
ing  plate.  For  an  end  slab  of  thickness  t,  the  value  of  v  at  the  side  in  question  for  boundary 
condition  (1.24)  is 

v  - -j3tanh(/3f)  (1.157) 

for  the  slab  just  below  the  upper  end  plate  and 

v  - +/3tanh(/3 1)  (1. 158) 

for  the  slab  just  above  the  lower  end  plate,  and  for  boundary  condition  (1.34)  it  is 


for  the  slab  just  below  the  upper  end  plate  and 


tanh(/3r) 


(1.160) 


for  the  slab  just  above  the  lower  end  plate.  A  practical  limitation  on  this  method  is  that  it  is 
unsafe  for  the  end  slabs  to  be  more  than  about  five  skin  depths  thick,  or  the  aforementioned 
numerical  instability  problem  can  render  the  use  of  (1.153)  dangerously  inaccurate.  Given 
these  end  slab  interface  values,  the  appropriate  matching  conditions,  and  expression  (1.150).  a 
complete  eigenfunction  may  be  constructed;  as  only  one  end  slab  interface  value  is  necessary 
for  this  purpose,  the  other  can  be  used  as  a  check  to  see  if  the  eigenvalue  has  been  properly 
determined. 

An  appropriately  sized  set  of  eigenvalues  can  be  calculated  to  a  good  first  approxima¬ 
tion  with  a  matrix  approach.  Equations  (1.20)  and  (1.29)  are  written  in  the  form 


U  /(.-)  -  kff(z) 


(1.161) 


L>  d 


. ±  _  v2(.) 

■y-’(z)  dz  (t 


lH  q(z )  -  kftq(z) 


(1.162) 


Lh--j t  -  yHz)  . 

dz ■ 

and  the  eigenfunctions  f(z)  and  q(z)  are  expanded  in  cosine  and  sine  series  respectively  so  that 
the  boundary  conditions  at  the  end  plates  are  automatically  satisfied; 


f(z)  -  a0  +  2L  cos(jnz/d) 

i- 1 


(1.163) 


q(z)  -  X.  bi  sin(Jnz/d) 


(1.164) 


where  d  is  the  distance  between  the  end  plates  and  for  the  sake  of  computational  practicality  the 
series  have  been  truncated  after  N  terms.  The  coefficients  are  then  represented  as  column  vec¬ 
tors  and  the  linear  operators  expressed  as  matrices; 


65 
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I  B„b,  =  kftb, 

i- 1 

A/,  —  ~  if  [Lt  cos{jn:/d))cos(ln:/ d)  dz  ,  0 

"  o 

^0/  =  -7  f  lit  cos(jwz/d)\  dz 

d  .1 

Bh  =  4  f  (£«  sin(jirz/d))sinUiTz/d )  dz 

a 


(1.166) 

(1.167) 

(1.168) 

(1.169) 


In  dealing  with  the  Z.(  operator,  one  should  note  that  the  second  term  of  the  operator  contri¬ 
butes  an  array  of  delta  functions  to  the  result  of  the  operation,  limit  calculations  indicate  that 
the  coefficient  of  this  term  is 


1  dyHz) 
y2(z)  dz 


Y,  8(z-r,)  log 


(1.170) 


where  there  are  N  homogeneous  layers  and  z,  is  the  coordinate  of  the  interface  on  top  of  layer 
number  i.  Due  to  the  truncation  of  the  trigonometric  series  representations  of  the  eigenfunc¬ 
tions,  the  series  cannot  properly  fit  the  most  sharply  curved  of  the  eigenfunctions,  and  as  a 
result  of  this  there  can  be  serious  errors  in  the  determination  of  their  associated  eigenvalues. 
The  only  effective  way  to  deal  with  this  problem  is  simply  to  bear  in  mind  that  some  of  the 
eigenvalues  computed  by  the  matrix  method  are  going  to  be  seriously  in  error,  with  those 
eigenvalues  associated  with  the  most  quickly  attenuating  modes  being  particularly  at  risk,  and 
be  prepared  to  detect  and  discard  the  bad  ones.  Square  matrices  with  dimensions  on  the  order 
of  80  are  easily  processed  for  their  eigenvalues  by  the  EISPACK  library  of  routines,  developed 
for  use  with  matrix  eigenvalue  problems  and  supported  by  many  major  scientific  computing 
facilities;  thanks  to  this  library  of  routines  it  is  quite  practical  to  use  outsized  series  expansions 
for  the  eigenfunctions  in  order  to  guarantee  that  a  useful  set  of  eigenvalue  solutions  remains 
after  the  grossly  erroneous  solutions  have  been  skimmed  off.  The  EISPACK  library  is  described 
in  detail  by  Smith  et  at.  (1976)  and  Garbow  et  al.  (1977)  in  the  two  primary  users'  guides  for 
this  library.  It  was  found  in  practice  that  a  useful  set  of  eigenvalues  could  usually  be  obtained 
by  having  the  number  of  terms  in  the  eigenfunction  expansion  to  be  double  the  number  of 
eigenvalues  desired  for  use,  thus  giving  double  the  desired  number  of  eigenvalues,  and  then 
discarding  the  eigenvalues  associated  with  the  most  rapidly  attenuating  modes  until  the  eigen¬ 
value  set  is  down  to  the  desired  size.  As  mentioned  previously,  the  rate  of  attenuation  of  a 
mode  is  roughly  proportional  to  the  absolute  value  of  the  imaginary  part  of  the  square  root  of 
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its  eigenvalue. 

As  well  as  rendering  some  of  the  eigenvalue  solutions  to  the  matrix  formulation  com¬ 
pletely  mendacious,  series  truncation  effects  can  easily  result  in  minor  errors  in  eigenvalue 
determination  that  are  nevertheless  sufficient  to  interfere  with  eigenfunction  construction  from 
the  eigenvalues;  a  method  is  needed  to  detect  the  serious  errors  and  refine  the  valid  if  approxi¬ 
mate  solutions.  The  method  that  was  used  for  this  purpose  is  to  construct  a  solution  for  v(z) 
from  the  top  layer  down  to  the  top  of  the  bottom  layer  and  determine  the  mismatch  between 
the  value  of  v  computed  at  this  point  from  above  and  the  value  computed  from  below,  this 
mismatch  being  negligible  in  magnitude  if  the  eigenvalue  used  for  the  construction  was  prop¬ 
erly  determined.  In  the  likely  event  that  the  mismatch  is  not  of  negligible  magnitude,  linear 
perturbation  methods  give  perturbation  of  the  mismatch  value  with  perturbation  of  the  pro¬ 
posed  eigenvalue,  and  a  version  of  the  Newton-Raiphson  method  adapted  for  use  in  the  com¬ 
plex  plane  may  then  be  used  to  home  in  on  the  proper  value  of  the  eigenvalue  in  question  pro¬ 
vided  that  it  is  not  too  far  away  in  the  complex  plane  from  the  approximate  value  calculated  by 
the  matrix  formulation.  In  the  event  that  the  trial  value  is  not  sufficiently  close  to  the  proper 
value  to  give  convergence  on  it.  the  procedure  will  still  try  to  find  a  valid  solution,  with  the 
most  probable  outcome  being  a  redundant  convergence  on  some  other  eigenvalue.  The  bottom 
of  the  topmost  layer  could  have  been  used  instead  of  the  top  of  the  bottommost  layer  for 
mismatch  calculation  and  reduction,  but  that  interface  will  typically  be  located  within  a  highly 
conducting  zone  modeling  seawater,  and  it  was  found  that  a  mismatch  parameter  calculated  at  a 
point  within  a  highly  conducting  region  can  lead  to  unusably  small  radii  of  convergence  in  the 
complex  k  plane  about  the  problem's  eigenvalues. 

In  the  case  of  a  receiver  located  at  the  interface  between  two  layers,  in  principle  the 
magnetic  field  and  the  horizontal  components  of  the  electric  field  should  be  calculated  as  having 
the  same  values  whether  the  receiver  position  is  taken  to  be  just  above  the  interface  or  just 
below  it;  in  practice,  the  method  of  Frieman  and  Kroil  sometimes  gives  somewhat  different 
field  values  for  the  two  receiver  positions,  the  cause  of  this  discrepancy  being  traceable  to  a 
numerical  problem.  If,  for  example,  the  interface  is  between  seawater  and  rock  with  an  electri¬ 
cal  conductivity  many  orders  of  magnitude  lower  than  seawater,  it  turns  out  that  the  calculation 
of  the  field  values  at  a  receiver  just  above  the  interface  involves  taking  the  differences  between 
large  numbers  to  get  a  result  smaller  by  several  orders  of  magnitude,  whereas  for  a  receiver  just 
beneath  the  interface  the  calculated  field  values  are  of  about  the  same  order  of  magnitude  as 
the  numbers  subtracted  from  each  other  to  get  them;  in  the  former  case  there  is  a  somewhat 
greater  opportunity  for  the  amplification  of  truncation  errors  than  is  present  in  the  latter  case, 
and  in  the  former  case  errors  in  the  magnitudes  of  the  calculated  electric  field  values  of  at  least 


its  eigenvalue. 


As  well  as  rendering  some  of  the  eigenvalue  solutions  to  the  matrix  formulation  com¬ 
pletely  mendacious,  series  truncation  effects  can  easily  result  in  minor  errors  in  eigenvalue 
determination  that  are  nevertheless  sufficient  to  interfere  with  eigenfunction  construction  from 
the  eigenvalues;  a  method  is  needed  to  detect  the  serious  errors  and  refine  the  valid  if  approxi¬ 
mate  solutions.  The  method  that  was  used  for  this  purpose  is  to  construct  a  solution  for  v(z) 
from  the  top  layer  down  to  the  top  of  the  bottom  layer  and  determine  the  mismatch  between 
the  value  of  v  computed  at  this  point  from  above  and  the  value  computed  from  below,  this 
mismatch  being  negligible  in  magnitude  if  the  eigenvalue  used  for  the  construction  was  prop¬ 
erly  determined.  In  the  likely  event  that  the  mismatch  is  not  of  negligible  magnitude,  linear 
perturbation  methods  give  perturbation  of  the  mismatch  value  with  perturbation  of  the  pro¬ 
posed  eigenvalue,  and  a  version  of  the  Newton-Ralphson  method  adapted  for  use  in  the  com¬ 
plex  plane  may  then  be  used  to  home  in  on  the  proper  value  of  the  eigenvalue  in  question  pro¬ 
vided  that  it  is  not  too  far  away  in  the  complex  plane  from  the  approximate  value  calculated  by 
the  matrix  formulation.  In  the  event  that  the  trial  value  is  not  sufficiently  close  to  the  proper 
value  to  give  convergence  on  it,  the  procedure  will  still  try  to  find  a  valid  solution,  with  the 
most  probable  outcome  being  a  redundant  convergence  on  some  other  eigenvalue.  The  bottom 
of  the  topmost  layer  could  have  been  used  instead  of  the  top  of  the  bottommost  layer  for 
mismatch  calculation  and  reduction,  but  that  interface  will  typically  be  located  within  a  highly 
conducting  zone  modeling  seawater,  and  it  was  found  that  a  mismatch  parameter  calculated  at  a 
point  within  a  highly  conducting  region  can  lead  to  unusably  small  radii  of  convergence  in  the 
complex  k  plane  about  the  problem's  eigenvalues. 

In  the  case  of  a  receiver  located  at  the  interface  between  two  layers,  in  principle  the 
magnetic  field  and  the  horizontal  components  of  the  electric  field  should  be  calculated  as  having 
the  same  values  whether  the  receiver  position  is  taken  to  be  just  above  the  interface  or  just 
below  it;  in  practice,  the  method  of  Frieman  and  Kroil  sometimes  gives  somewhat  different 
field  values  for  the  two  receiver  positions,  the  cause  of  this  discrepancy  being  traceable  to  a 
numerical  problem.  If,  for  example,  the  interface  is  between  seawater  and  rock  with  an  electri¬ 
cal  conductivity  many  orders  of  magnitude  lower  than  seawater,  it  turns  out  that  the  calculation 
of  the  field  values  at  a  receiver  just  above  the  interface  involves  taking  the  differences  between 
large  numbers  to  get  a  result  smaller  by  several  orders  of  magnitude,  whereas  for  a  receiver  just 
beneath  the  interface  the  calculated  field  values  are  of  about  the  same  order  of  magnitude  as 
the  numbers  subtracted  from  each  other  to  get  them;  in  the  former  case  there  is  a  somewhat 
greater  opportunity  for  the  amplification  of  truncation  errors  than  is  present  in  the  latter  case, 
and  in  the  former  case  errors  in  the  magnitudes  of  the  calculated  electric  field  values  of  at  least 
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50%  (with  a  precision  of  7  significant  digits  for  a  real  floating  point  number,  and  with  a  complex 
number  represented  by  two  real  floating  point  numbers)  were  not  uncommon.  Some  numerical 
exploration  of  the  problem  was  done  with  a  two  layer  model  with  seawater  on  top  and  poorly 
conducting  rock  below,  for  which  there  are  analytical  expressions  for  the  electric  and  magnetic 
fields,  and  agreement  between  these  expressions  and  the  method  of  Frieman  and  Kroll  was 
fairly  good  (usually  to  within  5%)  provided  that  the  receiver  position  used  with  the  latter  was 
below  the  interface.  It  was  noted  that  in  cases  where  the  conductivities  of  the  layers  above  and 
below  the  interface  differed  by  less  than  a  couple  of  orders  of  magnitude  the  discrepancies 
between  the  field  values  calculated  for  the  two  receiver  locations  were  usually  small  enough  to 
be  safely  disregarded. 
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Section  2 


Adaptation  of  the  Method  of  Frieman  and  Kroll 

to  Data  Inversion  Work 


Let  us  suppose  that  for  a  given  conductivity  profile  <rQ  we  have,  in  accord ’with  the 
method  of  Frieman  and  Kroll  (see  section  1),  determined  an  appropriate  set  of  eigenvalues, 
eigenfunctions,  and  Green’s  function  component  functions,  and  that  we  have  used  these  to  cal¬ 
culate  a  set  of  electromagnetic  field  values  at  selected  points  for  selected  frequencies.  In  the 
probable  event  that  we  find  these  field  values  calculated  from  the  profile  cr0  to  be  at  least  mildly 
inconsistent  with  experimental  measurements  of  the  electromagnetic  field,  we  wish  to  find 
some  conductivity  profile  cr\  which  yields  calculated  field  values  that  are  consistent  with  the 
experimental  measurements. 

The  first  step  in  the  inversion  algorithm  is  to  find  some  way  to  derive  expressions  of 

the  form 

U 

8  F,  =  X  g,(z)  8<r(z)  dz  ,  /  <■  1  to  N  ,  (2.1) 


where  F,  is  one  of  the  N  calculated  field  values,  6cr(z)  is  a  small  perturbation  of  the  conduc¬ 
tivity  profile  <t0(z),  8F,  is  the  accompanying  change  in  F,,  and  g,(z)  is  an  integration  kernal  to 
be  determined.  The  assumptions  of  vertical  stratification  and  horizontal  isotropy,  implicit  in  the 
method  of  Frieman  and  Kroll,  are  retained,  and  imposition  of  the  simplifying  assumption  that 
the  conductivity  profile  is  composed  of  discrete  homogeneous  slabs  turns  (2.1)  into 

SF,  -  £  C„  5<r,  ,  i  -  1  to  N  .  (2.2) 

/-i 


where  there  are  M  layers  in  the  profile.  Local  linearity  is  assumed  in  both  (2.1)  and  (2.2),  as 
global  linearity  between  <r  and  the  F,  does  not  exist. 

The  field  values  F,  are  calculated  from  a  collection  of  eigenvalues,  eigenfunctions,  and 
Green’s  functions,  which  in  turn  derive  from  the  solutions  of  the  equations  (1.20)  and  (1.29), 


d2f(z)  _  1  dy2(z)  df(z ) 

dz :  y2(z)  dz  dz 


L 


-  (y2(r)+*2)/U)  -0 


(2.3) 
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d2a(z) 


- 1 rv-->+**V-)  =o 


y2(r)  =  iu>iJ.Qcr(z)  (2.5) 

under  various  boundary  conditions.  Within  a  homogeneous  slab  both  (2.3)  and  (2.4)  have  the 


/"(:)- (y2+k2)f(z)  =  0  ;  (2.6) 

if  we  assume  f,  y 2,  and  k 2  to  be  the  perturbed  eigenfunction,  conductivity  profile  function,  and 
eigenvalue  respectively,  and  make  the  substitutions 

f(z)  -  f0(z)  +  8/0(r)  ,  (2.7) 

y2(z)  =  y<}(z)  +  SyJ(z)  , 

k2  «  k$  +  8 kfi  , 

where  /0,  -yj,  and  are  the  unperturbed  eigenfunction,  profile  function,  and  eigenvalue 
respectively,  we  get 

/o"  -  M+*J)/o<*>  +  8/0"(z)  -  (yS+kj  )8/0(z)  (2.8) 

-  (8y^+8^)/0(r)  -  (8y<j+8k<t)8f0(z)  -  0  . 

Assuming  that  the  perturbation  is  small  enough  so  that  it  is  safe  to  discard  the  second  order 
term,  we  then  have 

8fo"(z)  —  (y$+k$)8fo(z)  —  (8ytf+8ktf)f0(z)  ,  (2.9) 

since  in  order  for  f0(z)  to  be  a  proper  unperturbed  solution  we  must  have 

fo"(z)  -  (yM)fo(z)  -0  .  (2.10) 

The  general  solution  to  this  equation  is 

8/„(z)  -  A  cosh(Koz)  +  B  sinh (K0z)  (2.11) 

+  (8y$+8k$)  g(z)  , 

where  A  and  B  are  undetermined  constant  coefficients, 

Kt  -yt  +  k$  .  (2.12) 

and  g(z)  is  a  particular  solution  to  the  equation 


A 
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jf"(c)  -  K$  g(z)  -  f 0(r) 
Solutions  to  (2.10)  are  of  the  form 


(2.13) 


/,»(;)  -  C  cosh(Af0r)  +  D  sinMAfor)  , 


(2.14) 


whete  C  and  D  are  undetermined  constant  coefficients,  and  substitution  of  (2.14)  into  (2.13) 
gives  as  an  acceptable  solution  to  (2.13)  by  standard  mathematical  techniques  , 


g(:) 


C 


sinh(Koz) 

2/to 


+  D 


cosh(Koz)  -  —r  sinh(K0z) 
2A  o  2A< 5 


(2.15) 


Some  straightforward  if  tedious  algebra  gives  for  8/0(c)  from  (2.11) 

8/0(r)  -  E  8f0(z0)  +  F  8f0'(z0)  +  (8y«J+8/t02 )  G  (2.16) 

8/0'(r)  -  H  8/0(r0)  +  /  8/o'<r0)  +  (8-y<?-t-8 k$)  J  (2.17) 

E  -  cosh!*0<r— *«)}  (2.18) 

F  -  (l/A'o)  sinhj/f0(2-r0))  (2.19) 

G  -  g(z)  -  E  g(z0 )  -  F  g'(z0)  (2.20) 

H  -  Kq  sinh(tf0(z-r0)}  -  F  (2.21) 

/  «  cosh{A'0(z-zo)l  ”  E  (2.22) 

J  -  g'(z)  -  H  g(z0)  ~  1  g'(z0)  ,  (2.23) 


where  z o  is  within  the  same  slab  as  z.  These  relations,  with  the  appropriate  boundary  condi¬ 
tions,  can  be  used  to  construct  8/0(z)  and  solve  for  8kg  in  terms  of  the  conductivity  profile 
perturbations  So-,  for  the  M  homogeneous  layers  considered  to  make  up  the  profile.  It  is  help¬ 
ful  that,  although  special  methods  must  be  used  to  avoid  precision  problems  in  the  determina¬ 
tion  of  /o,  no  such  difficulties  arise  in  the  calculation  of  8/0  from  fo  once  /«  has  been  deter¬ 
mined. 

That  both  the  perturbed  and  unperturbed  eigenfunctions  must  satisfy  the  same  match¬ 
ing  conditions  at  the  layer  boundaries  allows  us  to  deduce  matching  conditions  for  the  perturba¬ 
tion  function.  From  the  boundary  conditions  (1.33).  (1.34),  and  (1.36)  we  have  that  perturba¬ 
tion  functions  on  solutions  to  (2.4)  must  satisfy  the  conditions  at  a  layer  interface  at  r-:() 
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8<7oi<-o)  =  8<7o2<-o)  (2.24) 

S^oi'(-o)  -  5<7o2'(-V  ,  <2.25) 

and  at  a  boundary  at  z=z0  between  a  layer  of  finite  conductivity  and  a  layer  of  finite  thickness 
and  infinite  conductivity 

8<?0(~o)  =  0  (2.26) 

From  the  boundary  conditions  (1.26)  and  (1.24)  respectively  we  have  that  perturbation  func¬ 
tions  on  solutions  to  equation  (2.3)  must  satisfy 


8/oi(-o)  “  8/02  (-o) 


(2.27) 


at  an  interface  at  z— z0  between  two  layers  of  finite  conductivity,  and 

5/o'(z0)  «  0  (2.28) 

for  a  boundary  at  2— z0  with  a  layer  of  finite  thickness  and  infinite  conductivity.  From  boundary 
condition  (1.23)  we  have  for  the  unperturbed  and  perturbed  functions 


/oi'(*o)  foi'(z o)  ,, 

“  W: 

/oi'(zq)  +  8/qi'(z0)  ^  /q2'(^q)  +  8/02' <zq)  (2  30) 

Wi  +  8y^,  y&  + 

using  the  binomial  expansion  and  discarding  terms  of  second  order  and  higher  gives  from 
(2.30) 


/o  1'  S/01  ’  /01 

2  '  2  +  4 

yoi  yoi  >oi 

_  /02'  +  8/02'  +  fp 2'SyJ; 

y$2  y&  >02 


and  use  of  (2.29)  then  gives 

Rr  ,  x?2  ,  y&  ,  , 

0/02  - r  ®7oi - r  7oi 

>01  >oi 


8y&  _  8yj| 
y<?2  ><li 


or,  alternatively. 


8/02’  “  ^~T  8/0/  —  /o2' 
>01 


8y&  _  8yj| 

>(?2  >021 


(2.31) 


(2.32) 


(2.33) 
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If  for  a  solution  to  (2.3)  wiili  its  boundary  and  matching  conditions  8/00  is  the  pertur¬ 
bation  function  value  at  the  interface  with  the  lower  infinitely  conducting  end  plate  (the 
assumption  of  infinitely  conducting  end  plates  capping  a  conductivity  profile  of  finite  thickness 
being  carried  over  from  section  1)  and  8/()W  is  the  value  at  the  interface  with  the  upper  end 
plate,  then  by  use  of  the  above  perturbation  function  continuation  and  boundary  conditions  we 
get  an  expression  of  the  form 

8/W  «  Co  8/oo  +  £  C,  tyl  (2.34) 

/-I 

+  Cw+i  8/fo  ”  0 

As  a  practical  measure  it  was  desirable  to  normalize  the  solutions  to  (2.4)  such  that 

</,/>  -  1  ,  (2.35) 

where  the  inner  product  has  the  definition 

<fufi>  - J  -  dz  ,  (2.36) 

and  requiring  (2.35)  to  hold  for  both  the  unperturbed  and  perturbed  eigenfunctions  gives  the 
condition  (again  discarding  terms  of  second  order) 


</o,8/0>  -  0  ;  (2.37) 

use  of  appropriate  continuation  and  boundary  conditions  and  some  straightforward  integration 
gives  from  this  an  expression  of  the  form 

</o’8/0>  -  D0  8/oo  +  £  D,  8yj,  (2.38) 

l-l 

+  Dw+i  8ko  “  0  , 

which  may  be  solved  simultaneously  with  (2.34)  to  yield  expressions  for  8 kj  and  8/oo  in  terms 
of  the  layer  conductivity  perturbations  alone;  for  example,  6k<$  may  be  expressed  in  the  form 

8 kj  -  £  K,  8 or,  ,  (2.39) 

i-i 

where  the  K,  are  complex  constant  coefficients.  These  in  turn  may  be  used  with  the  continua¬ 
tion  and  boundary  conditions  to  give  an  expression  for  8/0  of  the  form 

8/o(z)  -  ^  L,(z)  8(t, 


(2.40) 


The  determination  of  the  perturbation  functions  to  solutions  of  equation  (2.4)  with  its  associ¬ 
ated  boundary  conditions  is  handled  in  a  similar  manner. 

As  a  practical  matter  in  the  simultaneous  solution  of  equation  pairs  such  as  (2.34)  and 
(2.38),  it  was  found  that  a  straightforward  approach  of  scaling  followed  by  elimination  of  a  vari¬ 
able  through  subtraction  frequently  lead  to  precision  problems;  for  example,  often  one 
coefficient  would  be  many  orders  of  magnitude  greater  than  the  others  in  each  of  the  two  equa¬ 
tions,  and  it  would  be  the  equivalent  coefficient  in  both  equations.  In  practice,  it  was  found  that 
precision  problems  of  this  type  occurred  very  seldom,  if  ever,  in  work  with  perturbations  to 
solutions  to  equation  (2.4)  with  associated  boundary  conditions,  but  were  seldom  avoided  in 
work  with  equation  (2.3)  with  associated  boundary  conditions.  In  the  latter  case,  it  was  found 
that  if  one  multiplied  the  upper  end  plate  boundary  condition  equation,  of  which  equation 
(2.34)  is  an  example,  by  the  factor  F, 

- iliL -  (241) 

subtracted  the  result  from  the  norm  condition  equation,  of  which  (2.38)  is  an  example,  and 
used  the  result  in  the  simultaneous  solution  in  place  of  the  norm  equation,  precision  errors 
were  always  avoided.  The  functional  form  of  the  factor  F  was  determined  by  inspection  of  the 
coefficients  in  the  norm  equation  of  vanishing  terms  in  5 /0,w’.  and  the  numerical  factor  of  unity 
was  established  by  inspection  of  diagnostic  printouts  of  the  simultaneous  solution  routine. 

In  the  calculation  of  the  Green’s  function  component  functions  (see  for  example  equa¬ 
tions  (1.44)  and  (1.46))  it  was  found  desirable  as  a  practical  measure  to  normalize  each  com¬ 
ponent  function  so  that  the  function  value  was  equal  to  unity  at  the  layer  interface  nearest  the 
end  plate  boundary  from  which  the  function  was  developed,  this  normalization  condition 
together  with  the  appropriate  end  plate  boundary  condition  and  the  appropriate  value  of  k 2  was 
sufficient  to  uniquely  determine  each  component  function.  The  requirements  that  a  perturbed 
component  function  must  also  satisfy  the  condition  at  the  end  plate  boundary  and  be  equal  to 
unity  at  the  normalization  interface  are  sufficient  to  constrain  the  perturbation  function  to  the 
form 

8/c(z)  -  Co  8k2  +  £  C,(r)  So-,  ,  (2.42) 

1-1 

and,  since  the  appropriate  value  of  k 2  is  always  an  eigenvalue  of  some  eigenfunction  solution  of 
(2.3)  or  (2.4),  we  will  always  have  some  expression  of  the  form  (2.39)  to  substitute  into  (2.42) 
to  give  a  solution  of  the  form 
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8J-.(z)  =  £  D,(z)  8<r,  .  (2.43) 

<-i 

The  propagation  of  the  layer  conductivity  perturbations  from  this  point  on  to  the  field 
value  determinations  may  be  straightforwardly  accomplished  by  a  notational  convenience;  we 
define  a  vector-like  construction  such  that 

M 

A  =*  (Ao,A  i,^2'  *  *  “  Aq  4-  £  A,  b(T ,  ,  (2.44) 

i-t 

with  addition,  subtraction,  multiplication,  and  division  defined  such  that 

A  +  B  ”=  (A(\+Bq,A \+B\,A 2+S2,  •  •  ,Ay+ By)  (2.45) 

A  —  B  =“  (Aq—  Bq,A\— B\,A2~ B2,  ■  •  sA\f—By) 

A-B  =  (AqBq,AqB\+BqA  i,AoB2+BoA2,  •  • 

_y4  _  _  -4q^i  ^ 2  _  ^0^2  -4vr  _  AqB^ 

B  Bq  '  B0  B$  '  Bq  B<$  '  Bq 

The  multiplication  and  division  rules  are  derived  from  the  binomial  expansion,  with  terms  of 
second  order  or  higher  discarded;  note  that  it  is  a  property  of  this  notation  that  multiplication 
and  division  are  exact  inverse  operations  for  each  other,  as  they  are  for  normal  arithmetic. 
Replacement  in  the  field  value  computation  procedure  of  scalar  quantities  with  the  correspond¬ 
ing  vector-like  constructions  will  result  in  field  value  expressions  of  the  form 

F  -  (FlQ,F,hF, 2,  •  •  ■  ,F,„)  (2.49) 

-  F,0+  £  F,i  8o-;  . 

/-1 

The  electric  field  measurements  made  in  the  experiment  that  is  the  subject  of  this 
work  are  most  easily  represented  as  data  points  on  a  number  of  complex  planes;  unfortunately, 
due  to  experimental  difficulties,  absolute  phases  of  the  data  points  and  some  of  the  relative 
phases  were  for  practical  purposes  undetermined,  so  that  field  value  expressions  of  the  nature 
of  (2.49)  could  not  be  used  directly  for  data  inversion  work  on  the  available  data.  Instead,  an 
inversion  procedure  was  constructed  that  worked  with  the  lengths  of  vectors  connecting 
selected  pairs  of  data  points  in  their  complex  planes. 

r„  -  \F,-F,\  ,  (2.50) 

a  choice  of  parameterization  that  adequately  expresses  the  available  experimental  information 


(2.46) 

(2.47) 

(2.48) 
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while  being  unaffected  by  the  phase  uncertainty  problem.  The  notation  scheme  given  by  equa¬ 
tions  (2.44)  through  (2.48)  in  combination  with  the  standard  Pythagorean  formula  for  distance 
between  two  points  straightforwardly  gave  from  field  value  expressions  of  the  nature  of  (2.49)  a 
number  of  vector  length  perturbation  expressions  of  the  form 

/,  -  /,o  +  f  D„  8<r,  ,  (2. SI) 

i-i 

which  could  then  be  used  for  the  inversion  work. 

A  practical  difficulty  that  arose  at  this  point  was  the  fact  that  there  were  24  vectors 
whose  lengths  were  subject  to  fitting,  while  there  were  fewer  than  10  layers  in  the  conductivity 
profile  model  whose  conductivities  were  subject  to  perturbation  in  order  to  adjust  the  calculated 
vector  lengths;  a  least-squares  approach  was  indicated.  The  problem  was  expressed  in  the  form 
of  the  equation  system 

4  “  L  Ao  8<r>  ’  i  -  1  to  N  ,  (2.52) 

./“I 

where 

4-  w,  (/„-/, 0)  (2.53) 

and 

A,,  -  w,A,  .  (2.54) 

N  is  the  number  of  vector  lengths  to  be  fitted,  lie  is  the  experimental  value  for  the  i  th  vector 
length,  /l0  and  D„  are  as  in  equation  (2.51),  and  w,  is  a  weighting  factor  to  determine  the 
degree  of  emphasis  of  each  equation  in  the  least-squares  procedure.  The  weighting  factors  were 
chosen  to  be  proportional  to  the  standard  deviations  of  their  corresponding  vector  lengths,  the 
philosophy  being  that  in  the  fitting  procedure  those  vector  lengths  that  could  stand  to  have  the 
greatest  liberties  taken  with  them  are  those  that  are  the  most  poorly  determined  in  the  first 
place.  In  matrix  notation  we  can  express  the  equation  system  (2.52)  as 

rf-/4  6o-  ,  (2.55) 

where  d  is  a  column  vector  of  dimension  N,  So-  is  a  column  vector  of  dimension  M,  and  A  is 
anNXM  matrix.  The  least-squares  inversion  of  (2.55),  which  is  designed  to  find  the  So-  which 
minimizes  the  norm  of  the  difference  between  the  right  and  left  sides  of  (2.55).  employs  the 
method  of  singular  value  decomposition,  which  factors  the  matrix  A  into  the  form 
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A  =  U  Q  VT  ,  (2.56) 

where  U  is  an  N  X  M  matrix  such  that 

UTU  =  I  ,  (2.57) 

VT  is  an  M  X  M  matrix  such  that 

VT  V  =  V  VT  -  /  ,  (2.58) 

and  Q  is  an  M  X  M  matrix  all  of  whose  nondiagonal  elements  are  equal  to  zero.  We  then  have 
from  (2.55)  and  (2.56) 

8 v  V  R  UT  d  ,  (2.59) 

where  R  is  the  matrix  obtained  from  Q  by  replacing  the  nonzero  diagonal  elements  by  their 
reciprocals.  In  order  to  avoid  serious  injury  to  the  linearity  approximation  that  underlies  most 
of  the  work  leading  up  to  (2.52),  the  column  vector  d  is  scaled  down  to  such  a  size  that  none 
of  the  layers  of  the  conductivity  profile  model  have  their  conductivities  changed  by  more  than 
10%,  a  restriction  that  was  found  in  practice  to  be  sufficient  to  avoid  divergence  of  the  inver¬ 
sion  procedure  due  to  failure  of  the  linearity  approximation. 

A  useful  feature  of  the  method  of  singular  value  decomposition  is  that  it  can  be  made 
insensitive  to  machine  precision  error  and  other  sources  of  computational  noise.  In  a  problem 
such  as  is  represented  by  the  matrix  equation  (2.55),  there  will  often  be  one  or  more  conduc¬ 
tivity  perturbation  profiles  So-  that,  if  fed  through  A,  will  yield  column  vectors  d  of  small  or 
negligible  magnitude.  In  such  a  case,  when  the  problem  is  turned  around  and  a  data  vector  d  is 
input  to  yield  a  perturbation  profile  8 a,  noise  in  the  data  or,  in  extreme  cases,  truncation  error 
in  the  machine  that  processes  the  data,  can  result  in  one  or  more  of  these  null  perturbation 
profiles  being  grossly  exaggerated  (to  a  degree  roughly  proportional  to  the  negligibility  of  its 
associated  d  vector’s  magnitude)  and  superposed  on  the  desired  solution  profile;  this  noise 
amplification  effect  can  badly  obscure  the  physical  content  of  the  problem,  and  it  is  often  desir¬ 
able  to  completely  filter  such  null  perturbation  profiles  from  the  answer.  The  matrix  decomposi¬ 
tion  scheme  given  by  equations  (2.56)  through  (2.58)  points  out  the  null  profiles  buried  in  the 
matrix  A  and  makes  it  easy  to  eliminate  their  effects  on  the  solution  profile.  In  this  scheme.  U 
and  V  are  basically  coordinate  transform  matrices,  and  each  of  the  diagonal  elements  of  the 
matrix  Q  stands  for  one  of  a  number  of  mutually  orthogonal  perturbational  profiles;  if  such  an 
element  is  very  small  relative  to  the  diagonal  element  of  largest  magnitude,  the  profile  it 
represents  will  be  one  of  the  troublesome  profiles.  If  it  is  judged  that  this  profile's  contribution 
to  the  solution  profile  obscures  more  than  it  enlightens,  then  its  corresponding  element  in  Q 


can  be  set  to  zero  and  the  trouble  profile  will  be  deleted  from  the  basis  set  of  profiles  from 
which  the  solution  profile  is  constructed.  A  discussion  of  the  method  and  some  computer  rou¬ 
tines  that  implement  it  are  given  in  Wilkinson  and  Reinsch  (1971).  Although  the  need  for  the 
filtering  capability  of  this  method  never  arose  in  the  course  of  the  inversion  work,  it  was  reas¬ 
suring  to  have  this  capability  present  in  the  inversion  routines. 
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